home *** CD-ROM | disk | FTP | other *** search
/ Die Ultimative Software-P…i Collection 1996 & 1997 / Die Ultimative Software-Pakete CD-ROM fur Atari Collection 1996 & 1997.iso / g / gnu_c / pmlsrc23.zoo / pmlsrc / sin.c < prev    next >
Encoding:
C/C++ Source or Header  |  1994-03-19  |  6.9 KB  |  309 lines
  1. /************************************************************************
  2.  *                                    *
  3.  *                N O T I C E                *
  4.  *                                    *
  5.  *            Copyright Abandoned, 1987, Fred Fish        *
  6.  *                                    *
  7.  *    This previously copyrighted work has been placed into the    *
  8.  *    public domain by the author (Fred Fish) and may be freely used    *
  9.  *    for any purpose, private or commercial.  I would appreciate    *
  10.  *    it, as a courtesy, if this notice is left in all copies and    *
  11.  *    derivative works.  Thank you, and enjoy...            *
  12.  *                                    *
  13.  *    The author makes no warranty of any kind with respect to this    *
  14.  *    product and explicitly disclaims any implied warranties of    *
  15.  *    merchantability or fitness for any particular purpose.        *
  16.  *                                    *
  17.  ************************************************************************
  18.  */
  19.  
  20.  
  21. /*
  22.  *  FUNCTION
  23.  *
  24.  *    sin    double precision sine
  25.  *
  26.  *  KEY WORDS
  27.  *
  28.  *    sin
  29.  *    machine independent routines
  30.  *    trigonometric functions
  31.  *    math libraries
  32.  *
  33.  *  DESCRIPTION
  34.  *
  35.  *    Returns double precision sine of double precision
  36.  *    floating point argument.
  37.  *
  38.  *  USAGE
  39.  *
  40.  *    double sin (x)
  41.  *    double x;
  42.  *
  43.  *  REFERENCES
  44.  *
  45.  *    Computer Approximations, J.F. Hart et al, John Wiley & Sons,
  46.  *    1968, pp. 112-120.
  47.  *
  48.  *  RESTRICTIONS
  49.  *
  50.  *    The sin and cos routines are interactive in the sense that
  51.  *    in the process of reducing the argument to the range -PI/4
  52.  *    to PI/4, each may call the other.  Ultimately one or the
  53.  *    other uses a polynomial approximation on the reduced
  54.  *    argument.  The sin approximation has a maximum relative error
  55.  *    of 10**(-17.59) and the cos approximation has a maximum
  56.  *    relative error of 10**(-16.18).
  57.  *
  58.  *    These error bounds assume exact arithmetic
  59.  *    in the polynomial evaluation.  Additional rounding and
  60.  *    truncation errors may occur as the argument is reduced
  61.  *    to the range over which the polynomial approximation
  62.  *    is valid, and as the polynomial is evaluated using
  63.  *    finite-precision arithmetic.
  64.  *    
  65.  *  PROGRAMMER
  66.  *
  67.  *    Fred Fish
  68.  *
  69.  *  INTERNALS
  70.  *
  71.  *    Computes sin(x) from:
  72.  *
  73.  *      (1)    Reduce argument x to range -PI to PI.
  74.  *
  75.  *      (2)    If x > PI/2 then call sin recursively
  76.  *        using relation sin(x) = -sin(x - PI).
  77.  *
  78.  *      (3)    If x < -PI/2 then call sin recursively
  79.  *        using relation sin(x) = -sin(x + PI).
  80.  *
  81.  *      (4)    If x > PI/4 then call cos using
  82.  *        relation sin(x) = cos(PI/2 - x).
  83.  *
  84.  *      (5)    If x < -PI/4 then call cos using
  85.  *        relation sin(x) = -cos(PI/2 + x).
  86.  *
  87.  *      (6)    If x is small enough that polynomial
  88.  *        evaluation would cause underflow
  89.  *        then return x, since sin(x)
  90.  *        approaches x as x approaches zero.
  91.  *
  92.  *      (7)    By now x has been reduced to range
  93.  *        -PI/4 to PI/4 and the approximation
  94.  *        from HART pg. 118 can be used:
  95.  *
  96.  *        sin(x) = y * ( p(y) / q(y) )
  97.  *        Where:
  98.  *
  99.  *        y    =  x * (4/PI)
  100.  *
  101.  *        p(y) =  SUM [ Pj * (y**(2*j)) ]
  102.  *        over j = {0,1,2,3}
  103.  *
  104.  *        q(y) =  SUM [ Qj * (y**(2*j)) ]
  105.  *        over j = {0,1,2,3}
  106.  *
  107.  *        P0   =  0.206643433369958582409167054e+7
  108.  *        P1   =  -0.18160398797407332550219213e+6
  109.  *        P2   =  0.359993069496361883172836e+4
  110.  *        P3   =  -0.2010748329458861571949e+2
  111.  *        Q0   =  0.263106591026476989637710307e+7
  112.  *        Q1   =  0.3927024277464900030883986e+5
  113.  *        Q2   =  0.27811919481083844087953e+3
  114.  *        Q3   =  1.0000...
  115.  *        (coefficients from HART table #3063 pg 234)
  116.  *
  117.  *
  118.  *    **** NOTE ****      The range reduction relations used in
  119.  *    this routine depend on the final approximation being valid
  120.  *    over the negative argument range in addition to the positive
  121.  *    argument range.  The particular approximation chosen from
  122.  *    HART satisfies this requirement, although not explicitly
  123.  *    stated in the text.  This may not be true of other
  124.  *    approximations given in the reference.
  125.  *            
  126.  */
  127.  
  128. #if !defined (__M68881__) && !defined (sfp004)    /* mjr++    */
  129.  
  130. #include <stdio.h>
  131. #include <math.h>
  132. #include "pml.h"
  133.  
  134. static char funcname[] = "sin";
  135.  
  136. static double sin_pcoeffs[] = {
  137.     0.20664343336995858240e7,
  138.    -0.18160398797407332550e6,
  139.     0.35999306949636188317e4,
  140.    -0.20107483294588615719e2
  141. };
  142.  
  143. static double sin_qcoeffs[] = {
  144.     0.26310659102647698963e7,
  145.     0.39270242774649000308e5,
  146.     0.27811919481083844087e3,
  147.     1.0
  148. };
  149.  
  150.     struct exception xcpt;
  151.  
  152. double sin (x)
  153. double x;
  154. {
  155.     double y;
  156.     double yt2;
  157.  
  158.     if (x < -PI || x > PI) {
  159.     x = fmod (x, TWOPI);
  160.     if (x > PI) {
  161.         x = x - TWOPI;
  162.     } else if (x < -PI) {
  163.         x = x + TWOPI;
  164.     }
  165.     }
  166.     if (x > HALFPI) {
  167.     xcpt.retval = -(sin (x - PI));
  168.     } else if (x < -HALFPI) {
  169.     xcpt.retval = -(sin (x + PI));
  170.     } else if (x > FOURTHPI) {
  171.     xcpt.retval = cos (HALFPI - x);
  172.     } else if (x < -FOURTHPI) {
  173.     xcpt.retval = -(cos (HALFPI + x));
  174.     } else if (x < X6_UNDERFLOWS && x > -X6_UNDERFLOWS) {
  175.     xcpt.retval = x;
  176.     } else {
  177.     y = x / FOURTHPI;
  178.     yt2 = y * y;
  179.     xcpt.retval = y * (poly (3, sin_pcoeffs, yt2) / poly(3, sin_qcoeffs, yt2));
  180.     }
  181.     return (xcpt.retval);
  182. }
  183. #endif    /* __M68881__, sfp004    */
  184. #ifdef    sfp004
  185.  
  186. __asm("
  187.  
  188. comm =     -6
  189. resp =    -16
  190. zahl =      0
  191.  
  192. ");    /* end asm    */
  193.  
  194. #endif    sfp004
  195. #if defined (__M68881__) || defined (sfp004)
  196.  
  197.     __asm(".text; .even");
  198.  
  199. # ifdef    ERROR_CHECK
  200.  
  201.     __asm("
  202.  
  203. _Overflow:
  204.     .ascii \"OVERFLOW\\0\"
  205. _Domain:
  206.     .ascii \"DOMAIN\\0\"
  207. _Error_String:
  208.     .ascii \"sin: %s error\\n\\0\"
  209. .even
  210. | m.ritzert 7.12.1991
  211. | ritzert@dfg.dbp.de
  212. |
  213. |    /* NAN  = {7fffffff,ffffffff}        */
  214. |    /* +Inf = {7ff00000,00000000}        */
  215. |    /* -Inf = {fff00000,00000000}        */
  216. |    /* MAX_D= {7fee42d1,30773b76}        */
  217. |    /* MIN_D= {ffee42d1,30773b76}        */
  218.  
  219. .even
  220. double_max:
  221.     .long    0x7fee42d1
  222.     .long    0x30273b76
  223. double_min:
  224.     .long    0xffee42d1
  225.     .long    0x30273b76
  226. NaN:
  227.     .long    0x7fffffff
  228.     .long    0xffffffff
  229. p_Inf:
  230.     .long    0x7ff00000
  231.     .long    0x00000000
  232. m_Inf:
  233.     .long    0xfff00000
  234.     .long    0x00000000
  235.     ");    /* end asm    */
  236. # endif    ERROR_CHECK
  237.  
  238.     __asm(".even
  239. .globl _sin
  240. _sin:
  241.     ");    /* end asm    */
  242.  
  243. #endif    /* __M68881__ || sfp004    */
  244. #ifdef    __M68881__
  245.  
  246.     __asm("
  247.     fsind    a7@(4), fp0    | sin
  248.     fmoved    fp0,a7@-    | push result
  249.     moveml    a7@+,d0-d1    | return_value
  250.     ");    /* end asm    */
  251.  
  252. #endif    __M68881__
  253. #ifdef    sfp004
  254.     __asm("
  255.     lea    0xfffa50,a0
  256.     movew    #0x540e,a0@(comm)    | specify function
  257.     cmpiw    #0x8900,a0@(resp)    | check
  258.     movel    a7@(4),a0@        | load arg_hi
  259.     movel    a7@(8),a0@        | load arg_low
  260.     movew    #0x7400,a0@(comm)    | result to d0
  261.     .long    0x0c688900, 0xfff067f8    | wait
  262.     movel    a0@,d0
  263.     movel    a0@,d1
  264.     ");    /* end asm    */
  265.  
  266. #endif    sfp004
  267. #if defined (__M68881__) || defined (sfp004)
  268. # ifdef    ERROR_CHECK
  269.     __asm("
  270.     lea    double_max,a0    |
  271.     swap    d0        | exponent into lower word
  272.     cmpw    a0@(16),d0    | == NaN ?
  273.     beq    error_nan    |
  274.     swap    d0        | result ok,
  275.     rts            | restore d0
  276. ");
  277. #ifndef    __MSHORT__
  278. __asm("
  279. error_nan:
  280.     moveml    a0@(24),d0-d1    | result = +inf
  281.     moveml    d0-d1,a7@-
  282.     movel    #62,_errno    | NAN => errno = EDOM
  283.     pea    _Domain        | for printf
  284. ");
  285. #else    __MSHORT__
  286. __asm("
  287. error_nan:
  288.     moveml    a0@(24),d0-d1    | result = +inf
  289.     moveml    d0-d1,a7@-
  290.     movew    #62,_errno    | NAN => errno = EDOM
  291.     pea    _Domain        | for printf
  292. ");
  293. #endif    __MSHORT__
  294. __asm("
  295. error_exit:
  296.     pea    _Error_String    |
  297.     pea    __iob+52    |
  298.     jbsr    _fprintf    |
  299.     addl    #12,a7        |
  300.     moveml    a7@+,d0-d1
  301.     rts
  302.     ");
  303. # else    ERROR_CHECK
  304.  
  305. __asm("rts");
  306.  
  307. # endif    ERROR_CHECK
  308. #endif /* __M68881__ || sfp004    */
  309.